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An Updated Look at Binary Characteristics of Massive Stars in 

the Cygnus OB2 Association 

Daniel C Kiininki^^ & Henry A. Kobulnicky^ 

ABSTRACT 

O 

5^ ■ This work provides a statistical analysis of the massive star binary charac- 

teristics in the Cygnus 0B2 Association using radial velocity information of 114 
B3-03 primary stars and orbital properties for the 24 known binaries. We com- 
pare these data to a series of Monte Carlo simulations to infer the intrinsic binary 
fraction and distributions of mass ratios, periods, and eccentricities. We model 
^ ■ the distribution of mass ratio, log-period, and eccentricity as power-laws and find 

r^ ■ best fitting indices of a = 0.1±0.5, /3 = 0.2±0.4, and 7 = — 0.6±0.3 respectively. 

Y^- These distributions indicate a preference for massive companions, short periods, 

O , and low eccentricities. Our analysis indicates that the binary fraction of the 

^ ■ cluster is 44 ± 8% if all binary systems are (artificially) assumed to have P<1000 

days; if the power-law period distribution is extrapolated to 10'^ years, a plausible 
upper limit for bound systems, the binary fraction is ~90±10%. Of these binary 

K^ ■ (or higher order) systems, ~45% will have companions close enough to interact 

ly-^ ■ during pre- or post-main-sequence evolution (semi-major axis ^4.7 AU). The 

period distribution for P < 27 d is not well reproduced by any single power-law 
owing to an excess of systems with periods around 3-5 days (0.08-0.31 AU) and 
O ■ a relative shortage of systems with periods around 7-14 days (0.14-0.62 AU). 

We explore the idea that these longer-period systems evolved to produce the ob- 
served excess of short-period systems. The best fitting binary parameters imply 
that secondaries generate, on average, ~16% of the V-band light in young mas- 

;h . sive populations. This means that photometrically based distance measurements 

- - for young massive clusters & associations will be systematically low by ~8% 

(0.16 mag in the distance modulus) if the luminous contributions of unresolved 
secondaries are not taken into account. 

Subject headings: techniques: radial velocities — binaries: general — (stars:) 
binaries: spectroscopic — (stars:) binaries: {including multiple) close — stars: 
early-type — stars: kinematics — surveys 
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Introduction 



In clusters containing high concentrations of massive stars (B3 and earher), binary and 
multiple systems provide a virtual cornucopia of information about the massive stars and the 
nature of their formation process. In addition to the frequency at which they occur, we may 
examine the distribution of several quasi-preserved orbital parameters (e.g., period, mass ra- 
tio, and eccentricity) and gain insight into many open questions surrounding massive stars, 
such as whet her they form in isolation or as part of a competing menagerie of protostars (for 



a review see 



Zinnecker fc Yorkd 120071 ). Other issues include whether t here exists an evolu 



tionary link between massive bin ary systems and type Ib/c supernovae (JKobulnicky &: Fryer 
20071 : iPodsiadlowski et al.lll992l ). whether there is a correlation between the massive binary 
fracti on and the stellar density of the parent cluster (jSana et al.ll2008l : iGarcia &: Mermilliod 



20011 ). whether massive binaries are responsible for the runaway OB stars 



and what role these stars play in short and long 7-ray bursts f lFryer et al. 



Blaauwlll96ll ). 
1999h . Previ- 



ous and current endeavors to c haracterize the binary properties of massive stars include the 
galactic 0-star b inary study of iGarmany et al.l (Il980[ ). the Galactic 0- star speckle survey of 



Maiz et all (J2004J ). the multi-cluster 0-star binary statistical analysis of iGarcia fc Mermilliod 



(120011) . the VLT-FLAMES studies of the Large and Sm all Magellanic Clouds (lEvans et al. 



Sana et al. 



20061). the lucky imag in g 0-star s tudy of iMaizI (120101). the ada ptive optics 0-star studies 
of iTurner et"ai] J2OO8I ). IPuchenel (l200l[). and ISana et all (l201o[). the HST Fine Guidance 
Sensor survey of Cyg 0B2 by ICaballero-Nieves et al.l (I2OIII), and the in di vidual open cluste r 
binary spectroscopi c stud ie s of 



J2OO9I) . iMahv et al.l J2009I) 



Hillwid (l2006|),^_and 



Sana et al. (2008), 



Mermilliod] (1199511 . 



Rauw fc De Becked ( 120041 ). iDe Becker et al.l ( l2006h 
These studies have found a broad range of binary fractions (0-63%; ISana &: Evansll201ll ) and 
widely varying characterizations of the intrinsic orbital parameter distributions, including 
single power law descriptions with various indices and multiple-component descriptions re- 
quiring many free parameters. Reasons for the diverse, and seemingly incompatible, results 
in the literature stem from the varying methodologies and biases inherent in the surveys 
themselves. The majority of massive binary surveys in the literature include either a small 
sample of O stars within a given cluster or larger samples asser nbled from mult i ple smaller 



surveys encompassing different clusters fc associations (e.g., IGarmany et al.l Il980l : iGies 
19871 : lEvans et al.ll2006l : iMason et al.ll2009l ). The larger surveys have the advantage of better 
statistics, but the smaller surveys have the advantage of a more homogeneous sample (same 
formation environment and similar ages). 



The Cyg 0B2 radial v eloci ty survev (JKiminki et a. 



2OO9I . iKiminki et al.l I2OIOI . and iKobulnicky et al. 



2012 



.1 120071 iKiminki et al.ll2008l . iKiminki et al 



in prep; Papers I, II, III, IV, & VI), 
takes advantage of both a large and homogeneous sample of OB stars (114 massive stars) 
with radial velocity coverage of up to 12 years to probe the binary characteristics of massive 
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stars. iKobulnicky fc Fryerl (120071 KF07) analyzed the Cyg 0B2 radial velocity data from 
1999-2005 and compared the raw velocities with the expectations of Monte Carlo simula- 
tions over a range of mass ratio and orbital separation distributions and at varying binary 
fractionslll They find that th e likely binary fraction of OB stars within Cyg 0B2 is greater 
than 70% for an lOpikl (Il924j ) period distribution (i.e., log-flat) and mass ratio distribution 
best characterized by a power-law slope of a = —0.6 to 1.0 (depending on the choice of bi- 
nary fraction). Though KF07 are very th orough in addre ssing different possible distribution 
scenarios (a second "twin" population, a lHogeveenlll992l distribution. a lMiller fc Scaldll9'i 
secondary mass distribution, etc), they are forced to make certain simplifying assumptions in 
their code based on the limited binary information available. A few of the more significant 
assumptions are circular orbits, treating all higher order systems as binaries rather than 
triples, quadruples, etc, and finally, that all velocity variation stems from orbital dynamics 
rather than stellar atmosphere line profile variations, which may be present in massive stars. 
Assuming circular orbits can cause an underestimation of the binary fraction because fewer 
circular binaries need to be generated to simulate the seemingly singular eccentric systems 
that exist in the cluster. With regard to the next assumption, treating (difiicult-to-detect) 
triple and quadruple systems as binaries still yields a meaningful estimate of the binary 
fraction of the cluster, as long as this term is understood to include possible higher order 
systems. Lastly, the assumption that all velocity variation stems from orbital djTiamics can 
lead to an overestimation of the binary fraction or an overestimation of number of massive 
companions, but this is an unavoidable limitation to the KF07 analysis. 

We present here a follow-up to the initial results of Papers I, II, III, IV, & VI, and 
the earlier analysis of KF07. We use the results of nearly 12 years of spectroscopic obser- 
vations obtained on 114 massive stars in the Cyg 0B2 Association to model the observed 
mass ratio, orbital period, and eccentricity distributions composed from the orbital informa- 
tion of 22 close massive binaries in the cluster. We make use of all published close binary 
information for the cluster, both photometric and spectroscopic, including single-lined spec- 
troscopic binaries (SBls), double-lined spectroscopic binaries (SB2s), and eclipsing binaries. 
We also use the radia l veloc ity information of 110 OB stars from our original sample from 
Massey fc Thompson! (jl99ll ) to infer the binary fraction of the cluster. Additionally, we ex- 
amine, in more detail, an apparent exce ss of binary systems w ith periods between 3 and 5 days 
that have been o bserved in Cyg 0B2 ( iKiminki et al.l 120091 ) and other clusters (NGC 6231; 
Sana et al.ll2008l ). The current work provides an analysis using the largest sample of massive 



^We define the massive star binary fraction to mean tlie number of systems with two (or more) components 
divided by the total number systems containing at least one massive star. We caution the reader that other 
definitions exist in the literature. 
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stars and close massive binary systems ever compiled for a group of stars with a shared 
history. 

Section 2 reviews the data used in this work. Section 3 summarizes the method employed 
in this work to model the orbital period, eccentricity, and mass ratio distributions using a 
Monte Carlo approach. Section 4 presents the results of the Monte Carlo simulations and 
discusses the most probable intrinsic orbital parameter distributions and the probable binary 
fraction of massive stars within the cluster. Discussion in Section 5 explores possible physical 
origins for the observed orbital distributions and summarizes results that may inform models 
of close binaries as progenitors to supernova and 7-ray bursts. Finally, Section 6 summarizes 
the survey findings. 



The Data 



We obtained radial velocities for a total of 146 OB stars over 12 years. These 146 (given 
in Table 1 o f Paper I) were chosen from the UBV photometric and spectroscopic survey of 
Cyg 0B2 by lMassey &: ThompsonI (Il99ll ). Nearly half of the stars were previously identified 
as an OB-type and the other half had colors consistent w i th cla ssification of B3 or earlier. In 
this work, we utilize 110 of the 146 lMassey &: Thompson! (119911 ) stars, termed the "unbiased" 
sample because they were selected on a "binary-blind" basis (i.e., no previous evidence of 
a companion). The objects in the "unbiased" sample average 14 observations over 1999- 
2011 with a mean velocity uncertainty of 9.7 km s~^. Stars were removed from the original 
sample for having fewer than three observations at sufficiently high S/N needed to obtain 
reliable velocities (i.e., S/N<30 where velocity uncertainties are >40 km s~^) and for having 
spectral types later than B3. Additionally, we include the published orbital information for 
four sys t ems n ot in our original sample, but likely meml p ers of Cyg 0B2: Schulte 5 from 
Schultel Jl958h . A36, A45, and B17 from IComeron et al.l J2002h . This constitutes 114 OB 
stars and is designated the "complete" sample. The primary components in the "complete" 
sample have spectral types that range from B3V (e.g., MT298) to 03If (e.g., MT457) and 
masses that range from 7.6 M^r^ to 80 M^r^ (initial masses were interpolated from the stellar 
evolutionary models of iLejeune fc Schaereiil200ll ) with a mean mass of ~20 M© and mean 
temperature class of BO (71 B3-B0 stars and 43 09.5-03 stars). 

There are fewer stars than included in KF07 (who included stars from Table 5 of Paper I) 
owing to removal of seven stars having temperature classes later than B3 (MT196, MT239, 
MT271, MT453, MT493, MT576, and MT641), 13 stars having fewer than three observations 
with sufficient S/N (MT238, MT343, MT435, MT477, MT490, MT492, MT493, MT509, 
MT568, MT621, MT645, MT650, and MT712), and one star that is a newly turned-on Be 
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star (MT213). Additionally, ten stars not included in the analysis of KF07 but included 
here are Schulte 3, Schulte 73, MT005, MT179, MT186, MT222, MT241, MT267, MT421, 
and MT426. These ten stars were excluded from the earlier work because they had fewer 
than three observations prior to 2007. In total, this survey has produced full or partial 
orbital solutions for 19 binary systems. Combined with the other published data, there are 
24 known OB binary systems in Cyg 0B2. The locations of the first 20 are shown by the 
open circles in Figure 16 of Paper IV and listed in Table 6 therein. The remaining 4 are 
presented in Paper VI. 



3. Current Approach to Monte Carlo Analysis 

In keeping with prior works, we characterize the orbital parameter distributions as single 
power-laws. The mass ratio distribution is represented by f{q) oc q"" (where q = M2/M1) 
over a range of quiN to 1, the period distribution is represented by /(logP) oc (logP)^ over 
a range of periods from Pmin to Pmax days, and the eccentricity distribution is represented 
by /(e) oc e'^ over a range of cmin to cmax- In the current approach, we constrain the 
indices, a, /3, 7, by comparing the 12 years of radial velocity data, in addition to measured 
mass ratios, orbital periods, and eccentricities for the 22 binaries in Cyg 0B2 with periods 
under 30 days, to Monte Carlo-generated binary populations. For the standard application 
of the Monte Carlo (MC) c ode, we adopt qMiN = 0.0 05 (roughly representing the mass ratio 



of a B3V + M8V system; IPrilling fc Landoltl 120001 ). Pmin = 1 day, Pmax = 1000 days 



^MiN = 0.0001, and cmax = 0.9. The period upper limit is chosen based on our observing 
campaign time span of a few thousand days. The eccentricity upper limit is chosen based 
roughly on where an orbital system becomes unbound. We show later that our results are 
not sensitive to the exact choices for these upper and lower bounds, the main exception being 
that the inferred binary fraction does depend on the choice of Pmax- 

We created synthetic binary populations using Monte Carlo (MC) methods by randomly 
drawing periods, mass ratios, and eccentricities from power-law distributions with indices 
ranging from -2 to 2 in increments of 0.1, with binary fractions ranging from 20%-100% in 
increments of 5%. The MC code then generates a 4-dimensional parameter space composed of 
41 X 41 X 41 X 17 = 1, 171, 657 distinct "families" that are described by unique combinations of 
each parameter. During each iteration, radial velocities are calculated for a group of 110/114 
synthetic stars whose primary masses are determined from the adopted values for Cyg 0B2 
primaries. The number of synthetic velocities for each system is determined by the actual 
number of observations for each given star in the survey. The code samples the synthetic 
velocity curve at intervals corresponding to the actual observing cadence for that star. If the 
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system is, at random, designated a binary, a mass ratio, eccentricity, and period are drawn 
from the power-law distributions characterized by the family currently being examined. An 
inclination is randomly drawn from a probability distribution described by sin(i) oyer a n 



interval of to 7r/2, representing a random orientation on the unit sphere (iLaw et al.ll2009l ). 
The orbital phase and angle of periastron are randomly drawn from uniform distributions 
over their allowed values (i.e., 0-27r). The code adds Gaussian noise to each simulated radial 
velocity, where the magnitude of the noise is set by the observational uncertainty from the 
actual observations of the selected star. For example, when a system representing MT720 is 
simulated, a primary mass of 17.5 M0 is used and 32 radial velocities are calculated using 
the dates and errors listed in Table [TJ Finally, the code produces 300 realizations of each 
family, designated a "generation" . We compare the average properties of each generation to 
the Cyg 0B2 data. 

We constrain the binary fraction by comparing the x^ probability density function 
(PDF) for the simulated radial velocities to that of the unbiased sample of Cyg 0B2 stars. 
We use the same method described in Paper I, where the probability corresponds to the 
likelihood that the ~x^ for each star's radial velocity measurements, as considered around 
the weighted mean velocity, would be exceeded by chance given v = Nqbs ~ 1 degrees of 
freedom. We compare the observed x^ PDF to simulated ones for each generation using a 
2-sided Kolgorov-Smirnov (K-S) test. In our tests of the code, we were able to recover the 
binary fraction of a simulated dataset to within ±5%. 

In order to compare the Monte Carlo distributions of q, P, and e with the observed ones, 
it was necessary to apply a series of filters to the synthetic data to remove all binary systems 
with orbital solutions undetectable by this survey (i.e., the systems are designated "single" 
and the orbital parameters are not included in the final distributions). We first eliminate 
all synthetic systems have x^ probabilities >5%; such systems either have large errors, few 
measurements, and/or small velocity amplitudes that preclude full orbital solutions. The 
next filter removes all synthetic systems that have a semi-amplitudes Kmin <15 km s~^, the 
minimum semi-amplitude for which we can derive a complete orbital solution given typical 
velocity uncertainties. In exceptional cases where the stellar lines are narrow, we are able 
to do better. For example, MT174, with i^i = 9 km s~^ (Paper VI), has the smallest 
velocity semi-amplitude yet measured in the survey. Finally, we retain only systems that 
have periods within the 1.5-26.5 day range covered by the 19 systems with orbital solutions 
where our survey is highly complete (see Section 4.1). The final mean extracted synthetic 
orbital parameter distributions for each generation are then compared to the observed ones 
using 2-sided K-S tests. The final best-fit combination of a, /3, 7 is determined by product 
of the resultant individual K-S probabilities for each of the three parameter distributions. 
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In our testing of the code, we produced a series of simulated datasets (i.e., orbital 
parameter distributions and ephemerides) , similar to the Cyg 0B2 dataset, covering a large 
range of possible a, /3, and 7 combinations and binary fractions. The power-law indices 
ranged from —2.5 to 2.5 (and the MC code search grid was increased accordingly) and 
binary fractions ranged from 20% (the absolute minimum binary fraction of the Cyg 0B2 
sample based on current observations) to 100%. The simulated datasets were constructed 
using the observation cadences, velocity uncertainties, and primary masses from the Cyg 0B2 
data. The code was able to recover the original distribution indices to within ±0.3 rms and 
the original binary fraction to within 5% rms. The magnitude of the uncertainty on each 
index is limited by the current sample size of N=22 systems having orbital solutions. The 
uncertainty on the binary fraction is limited by the mean number of observations per object 
and the mean observational uncertainty. 



3.1. Limitations in the Code, Dataset, and Knowledge of Cyg OB2 

While our current approach explores a larger parameter space than KF07, incorporates 
eccentric orbits, and utilizes the information of 22 close massive binary systems in addition to 
the radial velocity information of 12 years of spectroscopic observations, there are still several 
limitations stemming from our modeling approach, our dataset, and our current knowledge 
of Cyg 0B2. 

Limitations owing to our modeling approach: As in KF07, we still consider all ve- 
locity variation as stemming from or bital dynamics rather than atmospheric phenomena 



such as pulsations. I Clark et al.l (120101 ) show that pulsations in evolved early-type stars can 
significantly affect radial velocity measurements, meaning that this effect may add scatter 
to radial velocity measurements in roughly one-quarter of our sample. We also still treat 
higher order systems as binaries. At present, we know of three triple or higher order systems, 
Cyg 0B2 N0.5, Cyg 0B2 N0.8, and MT429 (see references in Table 6 of Paper IV). In addi- 
tion, we limit characterization of each orbital parameter distribution to a single power-law, 
an assumption that may not hold. 

Limitations owing to the dataset: As mentioned, the Cyg 0B2 radial velocity survey 
probes binary systems with periods ranging from one day to a few thousand days, meaning we 
are not sensitive to extremely close (and extremely rare) systems with P <1 d. Therefore, 
the MC code only draws from period distributions with periods ranging 1 to 1000 days. 
Additionally, the vast majority of spectroscopic observations were obtained June through 
November when Cygnus is visible to northern latitudes. This, coupled with typical velocity 
uncertainties of 9.7 km s~^, reduces our sensitivity to long period, high eccentricity systems 



(i.e., systems that may pass through periastron and exhibit their largest velocity variation 
during an unobservable period), systems with extreme mass ratios, and systems with small 
inclination angles. Even though we have detected and computed solutions for systems with 
semi-amplitudes as small as 9 km s~^ (Paper VI), we attempt to circumvent these limitations 
by designating systems in the code with semi-amplitudes smaller than 15 km s~^ as "single". 

Another limitation comes in the form of not being able to detect unresolved SB2s. For 
systems exhibiting highly blended line profiles (even at quadrature), such as long period 
binaries with mass and luminosity ratios near unity, radial velocity variations of the blended 
line profile will be small and possibly undetectable at the r esolution of the sp ectra in this 
survey, thereby being designated as "single" by the code. ISana et al.l ( 120091 ) discuss this 
effect in their analysis of the massive star binary fraction of NGC 6611 and conclude that 
the effect is minimal (affecting their modeling of the binary fraction by only a few percent). 
Additionally, there is also a negligible bias toward detecting short-period, high-mass-ratio 
systems that show up as SB2s as it only takes one observation at quadrature to determine 
that the system warrants follow-up. These systems, however, have large velocity separations 
(i.e., large semi-amplitudes) and will get flagged as a "binary" by the code with only a few 
observations. 

Limitations owing to incomplete knowledge of Cyg OB2: A few studies have shown 
that Cyg 0B2 may be contaminated by an older population of stars, meaning that evo- 



lution ary effects may be present 



§,E^ 



Wright et al.ll2010l : lOrew et al.ll2008l : IComeron et al. 



20081 ). However, IComeron et al.l ( 120081 ) point out that the field star contamination is likely 



composed of F to late B stars at a radius of 1 degree. This exceeds the radial extent of the 
region examined in this study and concerns stars of a later type than are included here. On 
the other hand, some doub t remains in the mem bers hip of the 1 14 OB s tars included in this 
work as recent findings by IStroud et al.l (120101 ) and iLinder et al.l (120091 ) provide conflicting 
distance estimates to Cyg 0B2 based on their complete orbital solutions calculated for the 
eclipsing SB2s, Cyg 0B2 No. 5 (~1 kpc) and B17 (1.5-1.8 kpc) respectively. We maintain 
that the contamination is minor, if it exists, owing to the localized region we examine, but 
even the removal of a few periods or mass ratios could influence the distributions we observe 
at present. Similar to this is the assumption that the members of Cyg 0B2 have a common 
age. It has been suggested that eve n within the core reg ion of Cyg 0B2, there has been more 
than one epoch of star formation (jWright et al.ll2010l ). There are no clear predictions for 
how a population of close binaries should evolve over the 1-5 Myr lifetime of most massive 
stars in young clusters, so we have no way to model such effects or even know if they are 
important. 
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4. Results of the Monte Carlo Analysis 

We use the Monte Carlo code in conjunction with the "complete" sample of Cyg 0B2 
data to infer the underlying distribution of orbital parameters, that is, a, /3, 7. The orbital 
parameter information (but not ephemerides) for Schulte 9 and MT070 are excluded given 
their much longer periods. In order to address the binary fraction, however, we use the data 
of the "unbiased" sample. Ultimately, we find that the orbital parameter indices and binary 
fraction are essentially identical for both the "unbiased" and "complete" samples. 



4.1. Period Distribution 



Binary surveys that examine the orbital period distribution probe widely varying ranges 
in period/separation and stellar mass, and they find widely varying c haracterizations of its 
functi onal form. For solar-type stars with periods from 1-10,000+ days, iDuquennoy fc Mayor 
(ll99l[ ) find a log-normal distribution. Among massive star s, mo s t stud ies find a uniform 
distributi on in log space (0 = 0), sometimes known as an lOpikl (119241 ) distribution. For 
example, ICarmany et al.l (119801) find 0^0 for a large number of Galactic 0-star binaries in 



the northern hemisphere. iKouwenhoven et al.l (120071 ) also find an Opik's distribution for a 
survey of visual and spectroscopic bi naries in Sco 0B2 compos ed of A- and B-type primaries 
(for 0.7 days < -P < 3 x 10^ days). IZinnecker &: Yorkd (120071 ) speculate that rather than a 
simple power-law distribution for 0-star binaries, there may be a concentration of systems 
having periods be tween 2 an d 5 days. The 0-star binaries of NGC 6231 show this trend 
(ISana et al.ll2008l ). However, ISana &: Evans! (1201 ll ) conclude that a broken Opik distribu- 



tion provides a better description of the period distribution for spectroscopic 0-binaries in 
nearby clusters, where the broken distribution could better explain the higher concentration 
of systems with periods less than 10 days. 

The observed orbital period, eccentricity, and mass ratio distributions for the 22 known 
close massive bina ries in Cyg 0B2 are sh own in Figure [H The orbital period and mass ratio 
distributions from lCarmany et al.l (Il980[ ) are also shown and r epresented by the dott ed lines 
in the upper and lower panels, respectively. The Cyg 0B2 and lCarmany et al.l (Il980[ ) period 
distributions show a concentration between 3 and 5 days (with a peak around 4-5 days for 
Cyg 0B2). A possible explanation for such a concentration might be observational bias owing 
to the average observing run lasting a few days to a week. Three systems in Cyg 0B2 with 
periods between 3 and 5 days originate from other studies (MT421, MT429, and B17) where 
this may be the case. However, the remaining five systems in this period range originate 
from our study which has a time coverage of up to 12 years and individual observing run 
durations of 1-3 weeks. This allows for a sensitivity to orbital periods up to 25 days without 
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appreciable bias. This abundance of 3-5 day systems is not limited to the SB2s having the 
largest mass ratios either. Three of the nine SBls in Cyg 0B2 also have periods that are in 
the 4-5 days range. This is illustrated in Figure [21 which shows eccentricity versus orbital 
period (upper panel) and mass-ratio (lower panel) versus orbital period. The SBls and SB2s 
are the open and filled circles respectively, and the asterisks represent the orbital periods and 
eccentricities for the massive star binaries (primar ies of type B3 and earlier) in the Ninth 
Spectroscopic Binary Catalog ( jPourbaix et al.ll2004J ). One can easily see the concentration of 
points around 3-5 days and a peculiar shortage of systems with periods from 7-14 days in the 
Cyg 0B2 data. An initial hypothesis for this shortage of points might be that completeness 
plays a role. However, with multiple studies with varying observing cadences contributing 
to the catalog of binaries in Cyg 0B2, we maintain that a more likely observation would 
be to see a continually diminishing number of binary detections with increasing period as 
completeness declines. Given the existence of this excess in periods around 3-5 days, we 
caution that a simple power-law may be insufficient to describe the distribution of orbital 
periods for early-type stars in this cluster. 

In Figure [31 we show the results from the Monte Carlo analysis of the 22 close OB 
binaries in the "complete" sample. The four panels display the la (68.3%), 2(T (95.4%), and 
3cr (99.7%) contours relative to the maximum for the average probability of a generation as 
a function of two parameters for fixed best-fit values of the remaining two parameters. For 
example, the upper left panel shows the average probability as a function of a and 7 for 
fixed best-fit values of /3 (0.2) and B.F. (45%). The upper right panel similarly shows the 
average probability as a function of /3 and 7 for fixed best-fit values of a (-0.1) and B.F. 
(45%). The cross in each panel marks the location of the best fit for each parameter (i.e., the 
maximum of the average probability). We find that the period distribution is confined at the 
la level to —0.3 < /3 < 0.5, with a best fit of /3 = 0.2 ± 0.4. This points to a nearly log-fiat 



distributio n, consistent with Opik ' s Law (IOpiklll924l). c onsist e nt with the pre v ious -star 
findings of iGarmany et al.l (Il980l ) , iKouwenhoven et al.l (120071 ) , ISana fc Evans! (1201 ll ) , and 
KF07. Figure [H shows the cumulative period distribution of the Cyg 0B2 data (diamonds) 
and the best fit cumulative period distribution from the MC code (solid curve). The hump 
in the distribution at log{P) ~ 0.65 or P ~ 4.5 days greatly infiuences the fit shown, and 
it is clear that a power-law of /3 = might be a better fit if this hump were not present or 
less pronounced. Indeed, no single power-law provides a good fit to the cumulative period 
distribution. We address the peculiar nature of the period distribution in Section 5. 

Figure [D shows the survey completeness as a function of orbital period for the best- 
fitting power-law indices. We use the MC code to estimate the survey completeness by 
tallying the fraction of synthetic systems deemed detectable in our survey given the time 
sampling and velocity semi-amplitude limit of Kmin = 15 km s~^. The survey is 88% 
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complete at P = 10 days and > 83% complete at P = 30 days. Our completeness analysis 
implies that, given 14 detected binaries with periods P < 10 days, there are an additional 
1-2 undetected binary systems in this range. Statistically, the Monte Carlo code shows that 
all of these would be undetectable by the survey because they have some combination of 
low inclination, small mass ratio, large eccentricity, and larger radial velocity uncertainties. 
If, for the sake of argument, we ignore the excess of short-period systems and assume that 
the /3 ~ power-law holds, the ~15 binaries in the 1-10 day period range would imply 
an equivalent number between 10 and 100 days. There are 8 known systems in this range, 
leaving an additional 7 putative binaries as undetected. The mean completeness over this 
range is ~ 82%, suggesting that ~12 of the 15 would be detectable by our survey. This 
implies that the ongoing Cyg 0B2 survey should uncover as many as ~5 additional binaries 
with periods between 10 and 100 days, given sufficient long-term data. For periods longer 
than 100 days, w here there are presently only two systems with orbital solutions, Schulte 9 
(JNaze et al.l 120101 ) and MT070 (Paper VI), detection becomes more difficult. In the case of 
Schulte 9, with a period of ~2.355 years and an eccentricity of nearly 0.7, even an observing 
campaign such as this one experienced difficulty in uncovering this massive "twin" system. 
All of the foregoing discussion, however, assumes an Opik period distribution, and we have 
already shown in Figure H] that no single power-law adequately reproduces the data at the 
shortest periods, P < 14 days. 



4.2. Eccentricity Distribution 



Observationally, binaries with solar-type primaries and periods lon ger than a year have a 
broad distribution of orbital eccentricities with a median value of 0.55 (JDuquennoy fc Mayor 
19911 ). For periods shorter than this, these systems have a higher concentration near zero 
eccentricity stemming from circularization of the orbits. This is also shown to be true 
for 0-star binaries with periods as long as 10,000 days (ISana fc Evand l201ll ). Although 



the method by which binary orbits are circularized and the critical period below which 
all are circular is unclear from existing data, a standard eccentricity versus period scatter 
diagram, like Figure |2l illustrates that binaries with short periods do tend toward nearly 
circular orbits . The figure shows a decline of mean eccentricity for periods under 10 days. 
Giuricin et al.l ( 1l984j ) found that nearly all (92%) early-type binaries with periods shorter 
than two days had an eccentricity less than 0.05. They also found that the ec centric i ties o f 
nearly all of the early-type binaries in their survey could be explained by the IZahnI (119771 ) 
tidal theory, which describes the circularization time for an equal mass binary with an 
initial eccentricity of e < 0.3 a s being related to the main sequenc e lifet i me of that binary. 
KhaliuUin fc KhaliuUinal ( 120101 ) provide a modified version of the IZahnl (J1977I ) theory and 
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address the evolution of the orbital parameters over the circularization time of the orbit . 
The circularization timescales are, on average, a third of tho se predicted by IZahru (jl977| ) 
and in bett er agreement than the predictions of lTassoull (Il988l ) with the massive star binary 
statistics of iGiuricin et al.l (Il984f ). The predicted circularization time for an early-type (pair 
of BOV stars), nearly equal- mass binary at a separation comparable to the Cyg 0B2 early 
B-type binaries (e.g., asini ~0.16 AU for MT720) is c lose t o the main-sequence lifetime of an 



early B star. However, a study by lWitte fc Savonijd ( 1200 ll ) found that there are a significant 
number of early- type binaries with periods up to 10 days with an eccentricity near zero that 
cannot b e explained with the IZahnI (119771) theory and are unlikely to have been primordially 
circular. iKrumholz fc Thompson! ( 120071 ) have proposed a method by which orbits become 
significantly circularized during the star formation process, thereby diminishing the conflict 
between circularization time and the theoretical life span of massive stars. However, the 
theory assumes the binary is formed concurrently with the formation of the components, 
and that is a highly debated topic in massive star formation theory ( IZinnecker fc Yorke 
2007h . 



In Figure E] we show the MC code best fit to the Cyg 0B2 eccentricity distribution. 
As in Figure SJ diamonds show the normalized cumulative distribution (left Y-axis) and the 
solid curve shows the normalized cumulative distribution for 7 = —0.6. Figure |3] (upper left 
and right panels) shows that the index 7 is confined, at the la level, to —1.0 < 7 < —0.4 
with a best fit of 7 = —0.6 ± 0.3. The steep negative power-law is required to reproduce the 
abundance of systems with low eccentricity. A K-S test shows that the 7 = —0.6 power law 
is compatible (i.e., the observed sample and the synthetic sample are consistent with being 
drawn from the same parent distribution) with the data at the ~90% level. The absence of 
eccentricities with e > 0.55 i s similar to what is observed for short-period binaries in other 
OB clusters (e.g. NGC 6231: ISana et al.ll2008l ). Highly eccentric systems are apparently rare 
in this period regime, either because circularization mechanisms are efficient, or possibly 
because they are more easily disrupted during encounters with other cluster stars. Finally, 
we note that the eccentricity distribution within Cyg 0B2 is also compatible with that of 
OB binaries in the 9th Spectroscopic Binary Catalog at the ~50% level. 

We show an estimate of the survey completeness as a function of eccentricity for two 
different period ranges in Figure [3 The solid and dot-dash curves represent the completeness 
for P < 1000 days and P < 30 days respectively. The plot indicates that the survey is ~75% 
complete at e ~ 0.65 for binaries with P < 1000 days and ~88% complete at e ~ 0.75 for 
binaries with P < 30 days. For P < 30 days, the survey is likely to have detected the vast 
majority of systems, even those with large eccentricities, if they exist. 
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4.3. Mass Ratio Distribution 



The nature of the mass ratio distribution is one of the more highly debated and un- 
resolved questions pertaining to massive binaries and their formation. The distribution is 
difficult to characterize because of the unavoidable bias of observational studies, including 
this one, against low-g systems. The general consensus is that short-perio d binaries have 



a diffe rent mass ratio distribution corn pared to the long-period binaries. iGoldberg et al. 



(1199 ll ) and iDuquennoy fc Mayon (119911 ) show this with lower-mass binaries, but they also 
point out that the period range over which this change occurs is still ill-defined. Most stud- 
ies on the mass ratio distri bution are for mid to late-type stars. A few examples include 
Duquennoy fc Mayorl (ll99ll ). who find a < — 2 (for q > 0.1 and periods 1-10,0 00-1- days) 
for a large survey of solar- type spectroscopic binaries, iKouwenhoven et al.l ( 120071 ) . who find 
a = —0.4 to —0.5 for a large survey of spectroscopic and visual A- and B-type primaries in 
Sco 0B2, IPisher et al.l J2005I ). who find a ~ for all th e SBls and a > for all the SB2s 
in the solar neighborhood (for periods 1-500 days), and iHogeveenl ( 119921 ). who find a ~ 
for mass ratios of g < 0.3 and a < —2 for mass ratios of g > 0.3 in binaries of B-type 
an d later (for periods <1 day t o 10,000-1- days). Among the massive stars, studies such 



as 



Pinsonneault fc Stanekl (120061 ) conclude from their study of 21 eclipsing massive binaries 



(P < 5 days) in the Smal l Magella r iic Cl oud, that the secondary mass function for massive 
binaries does not follow a ISalpeten ( 1l955l ) slope. Instead, the authors find a fiat distribution 
of mass ratios and predict that rough ly 55% of all massive systems are part of a "twin" 
system with a mass ratio of g > 0.95. ISana &: EvansI ( 120111 ) have refuted this proposition, 
but also find a preference for more massive companions among 0-stars in clusters. 

The distribution of observed mass ratios in Cyg 0B2, shown in Figure (H appears rela- 
tively fiat, except for a slight preference for mass ratios between 0.7 and 0.9. This may be an 
artifact of using a standard histogram for small numbers. Not surprisingly, the upper end of 
the distribution is populated exclusively by SB2s, while the lower end of the distribution is 
almost entirely populated by SBls. For five of these single-hned systems, we calculated the 
expected orbital inclination, (i), between imin and 90°, using the inclination lower limits 
estimated from Paper II, Paper III, and Paper VI. We calculate the most probable mass ra- 
tios to be 0.27, 0.30, 0.03, 0.22, and 0.11 for MT059, MT145, MT174, MT258, and MT267, 
respectively. These are close to the lower-limit values for q because (i) is not far from the 
upper limit on i. There is a possibility, of course, that for one or more of these systems, the 
orbital inclination could be less than (i), and therefore, the system(s) would have a more 
massive companion. In such a case, our present estimates may incorrectly infiate the number 
of systems with mass ratios between 0.05 and 0.30, and the true distribution would not be 
as fiat as currently shown in Figure [TJ 
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From the MC code we find a best fit of a = 0.1 ± 0.5 for tlie mass ratio distribution. 
Figure [8] sliows the normalized cumulative distribution of mass ratios (diamonds), along 
with the corresponding simulated cumulative distribution for a = 0.1 (thick curve). The la 
contours (top left and bottom left panels in Figure E]) show a range of —0.3 < a < 0.7. Such 
a shallow positive power-law indicates that mass ratios in the cluster are quasi-uniformly 
distributed between q = 0.005 and q = 1.0 with a preference for la rger mass ratio s. This 



type of distribution implies that the secondaries are not drawn from a ISalpeted (119551 ) stellar 
IMF, in contrast to lower-mass binary populations which exhibit an abundance of low mass 
ratios. Furthermore, the histogram of m ass ratios in Figure El i s also inconsistent with 



the "twin-heavy" population suggested by iPinsonneault fc Stanekl (120061 ). whereby 45% of 
systems are predicated to have q > 0.95. 

Figure M shows the completeness of the survey as a function of mass ratio. For systems 
with P < 30 days (dot-dash curve), the survey is at least 90% complete down to g ~ 0.25, 
but there is a sharp drop in the fraction of systems detected for q < 0.25. The average 
completeness is 10-20% lower, overall, for longer period systems with P < 1000 days (solid 
curve). Our measurement of MT174 (g = 0.03; Ki = 9 km s~^) below our nominal semi- 
amplitude detection threshold of 15 km s~^ provides evidence that such low-amphtude, low-g 
systems exist and are detectable in our survey. Figure |9] shows that incompleteness is an 
issue only at the lowest mass ratios. Any attempt to correct for incompleteness at g < 0.25 
would add low-g systems to the histogram and serve to fiatten the distribution. 



4.4. Binary Fraction 

Table [2] shows a selection of the larger surveys that investigated the 0-star binary 
fraction within individual clusters/associations. Column 1 contains the common cluster 
identification, column 2 is the number of O stars included in the survey, column 3 is the 
number of early B stars included in the survey, column 4 is the lower limit on the 0-star 
binary fraction, column 5 is the lower limit on the massive star binary fraction (O and early-B 
primaries), and column 6 references the original works. The observed 0-star binary fractions 
(designated b.f. rather than B.F. here to distinguish from the intrinsic binary fraction) within 
the clusters listed range widely from - 67% (column 4). More inclusively, the massive star 
b.f. (with primaries B3 and earlier; column 5), varies between 4% (NGC 330 in the SMC) 
and 53% (NGC 6231). In Paper I, we find an initial lower limit on the binary fraction of 30% 
based on the x^ probability density function for radial velocities from the first six years of 
observations of Cyg 0B2, which would place it near the middle of the binary fraction range 
among massive stars. KF07 find a more probable binary fraction of B.F. > 70% based on 
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their Monte Carlo analysis of the raw radial velocities, which would place it closer to that of 
NGC 6231. However, based purely on the binary detections of Papers II-IV & VI, the hard 
lower limit on the binary fraction among massive stars in Cyg 0B2 is 21% (i.e., 24/114). 

We show the observed x^ probability density function (PDF) for the radial velocities of 
Cyg 0B2 stars in Figure [TOl There are 53 having P(x^) '^) <5% and 43 with Pix^, v) <1%. 
These systems encompass the 24 known binary systems and additional candidates for radial 
velocity variables. The number of variables has increased from Paper I despite removing 
some low S/N spectra that may have caused extraneous velocities. This is understandable 
given the additional data acquired on the sample, especially on long period binaries like 
MT070 (P = 7.58 years). Figure [U] shows the MC code best fit to the cumulative x^ PDF 
for radial velocities in Cyg 0B2 using the "unbiased" sample of 110 stars. The solid curve 
represents the cumulative x^ PDF produced by the MC code with best-fit indices of a = 0.1 
(gjv/77v=0.005), [3 = 0.4 (Pa/^x=1000 days), and 7 = -0.6 (eA/w=0.0001, eMAx=0.9). The 
open diamonds represent the normalized cumulative x^ PDF for the unbiased sample of OB 
primaries. The best MC code fit results in a binary fraction of B.F. = 45 ± 8% if all binary 
systems are assumed to have P < 1000 days. Simulated binary populations with binary 
fractions of 40-50% in increments of 1% and with fixed a, /3, and 7, were individually 
examined for better agreement with the data. These additional simulations show that a 
binary fraction of 44% yields a slightly better fit to the data. 

Figure [T2] shows the results of the Monte Carlo analysis for the 18 OB close binaries 
comprising the "unbiased" sample. Similar to Figure [31 the four panels provide the la, 2a, 
and 3a probability contours for a, /3, 7, and binary fraction. The three panels depicting 
probabilities for power-law indices (upper left, upper right, lower left) are nearly identical 
to the results from the MC analysis of the "complete" sample shown in Figure [31 The 
best-fit P is slightly larger (0.4 instead of 0.2) owing to the exclusion of several shorter- 
period binaries from the unbiased sample. The best-fit a and 7 are the same, just with 
larger uncertainties. The contours in the lower right panel show that the binary fraction is 
confined to 37% < B.F. < 49% at the la level and 30% < B.F. < 60% at the 3a level for 
Pmax = 1000 days. This result is consistent with the range of massive binary fractions for 
Galactic clusters listed in Table [H In order to compare our current results to the Monte 
Carlo analysis of KF07, who explored a larger range of periods (separations of 0.02-10,000 
A.U. corresponding to P ~1-10^ days for two 20 M© stars), we also performed a MC run 
with Pmax = 10^ days. In both analyses 10,000 A.U. is considered a reasonable upper hmit 
for bound systems. In this case, we find the binary fraction must be ~ 100% in order to 
replicate our observed x"^ PDF. This is in rough agreement with the conclusions of KF07, but 
based on extrapolating the power-law distributions of mass ratios and periods determined 
from close binaries over many orders of magnitude to characterize loosely bound systems. 
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Is an extrapolation of the power law index /3 to periods as long as 10^ days warranted 
when our radial velocity survey is se nsitive only to systems with periods less than about 10 
years? ICaballero-Nieves et al.l ( 120111 ) conducted an adaptive optics and HST Fine Guidance 
Sensor survey of 75 Cyg 0B2 systems and detect highly probable companions around 31 OB 
primaries at separations >200 A.U. (P>630 years). Of these, nine (9/75 = 12%) have 
probable periods in the range 10^-10^ years (S. Caballero-Nieves, private communicationjj. 
If we adopt Pmax = 10^ years (~ 365,250 days), the MC code provides a best fit binary 
fraction of ~90%, and this extrapolation of the /3 ~ 0.2 power law predicts approximately 
15 binaries with orbital periods in each 1 dex interval (i.e., 15 with periods between 10 yr 
and 100 yr, 15 between 100 yr and 1000 yr, etc.) given the 110 stars in our unbiased sample. 
This translates into about 14% of the binary systems falling within each logarithmic period 
interval. Hence, the predictions of an extrapolated power-law appear consistent with the 
emerging high angular resolution observations at larger separations. We therefore suggest 
that the binary fraction among massive stars may be as high as 90 ± 10% if periods as long 
as 10"^ years are allowed. 



4.5. Sensitivity to other model assumptions 

In addition to exploring the effect of varying Pmax, we also examined the effect of 
varying Qmin and the semi- amplitude detection threshold, Kmin, in the MC code. As Kmin 
is changed from 15 km s~^ to 20 km s~^ or larger, a greater fraction of synthetic systems are 
deemed undetectable by the code, and the detected systems have smaller separations and/or 
larger companions. Therefore, we expect some combination of increase in a (larger mass 
ratios), decrease in /3 (shorter periods), and increase in binary fraction. We find insignificant 
change (i.e., a, /3, 7 are the same) when Kmin = 20 km s~^. A significant change in 
the indices is not observed until Kmin is increased to much larger values of ~ 50 km s~^. 
We are also confident our detec tion threshold is rnuch closer to 15 km s~^, as evidenced 
by the measurement of MT070 (JKobulnicky et al.l 120121 ). We conclude that the 



semi-amplitude for binary detection is appropriately set at 15 km s 



minimum 

1 



In the case of increasing Qmin, the lowest mass ratio simulated in the code, we primarily 
expect to see a decrease in a and a decrease in the binary fraction. Binaries will have more 
massive companions and semi-amplitudes will be larger, meaning that a smaller fraction of 
binary systems will be deemed undetectable in the code. We observed only a small decrease 



^We use a distance of 1700 pc, correct statistically for projection effects, and assume a total system mass 
of 30 M0 to turn the angular separations into estimated periods. 
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in binary fraction in runs with Qmin = 0.1 (a few percent and still consistent with our initial 
results). We also find only a small increase in a and no change in /3 or 7. However, increases 
in a become larger when Qmin > 0.1. For example, a changes from 0.1 to —0.4 for the 
"complete" sample when Qmin = 0.2; however, such a large Qmin is probably unphysical 
given the existence of MT070 (g > 0.03) and MT267 {q > 0.11). These results, for both 
tests, indicate that our MC results are not sensitive to small changes in the adopted values 
of the mass ratio lower limit or semi- amplitude threshold. 



5. Discussion 

5.1. Possible Origin of the Short-Period Excess 

The presence of an excess of binary systems with periods around 3-5 days (with a peak 
at 4-5 day s) and a sho r tage o f systems with periods around 7-14 days is not only observed in 



Cyg 0B2. ISana et al.l ( 120081 ) observe a similar occurrence in NGC 6231, and the excess and 
shortage of systems occur at roughly the same ranges in period. Systems in Cyg 0B2 that 
lie within the excess encompass a range of spectral types, mass ratios, and binary types (i.e., 
eclipsing, SBl and SB2) and are scattered randomly across the cluster, having no preferred 
location either in dense or sparse regions. We do not believe this peculiar distribution can 
be attributed to observational biases given the extended time coverage of this survey and 
the different observing cadences of the observing runs that have contributed to the Cyg 0B2 
binary catalog. While the existence of an excess of systems with periods of around 3-5 days 
seems well-established in the literature and in our present data, the possibility of a shortage 
in the 7-14 day range has not been previously no ted. A comparison be tween the Cyg 0B2 



cumulative period distribution with the data from lGarmany et al.l ( 1l980l ) reveals a flattening 
in this 7-14 day range (though we should be cautious given the use of a standard histogram 
and small numbers). Their survey (of 28 systems) contains three in this 7-14 day range 
while ours (22 systems) contains one. 

While we find that no single power-law can adequately describe the log-period distri- 
bution in Cyg 0B2, it is unlikely, given any power-law distribution, to find so few systems 
between 7 and 14 days. We examined the probability that only a single system would 
be observed in the 7-14 day period range given a random drawing of 22 systems with 
1.46 < P < 25.2 days (for /3 =-1 to 1). Random samples were drawn 10,000 times and 
the final percentages were adjusted for a 90% completeness level (the mean level as mea- 
sured in Figure for the 7-14 day range). The likelihood of detecting a single system for 
/3 = is ~2%. However, if the power-law describing the log-period distribution is assumed 
to be closer to —1, then the probability rises quickly to ~13% (with ~ 5% for /3 = —0.5, and 
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< 1% for P =0.5 and 1). This tells us that if the power-law(s) that describe the log-period 
distribution have /3 > —0.5, the shortage is a real phenomenon at the 2a level. 

In an attempt to replicate the Cyg 0B2 period distribution with a physically motivated 
prescription, we made a version of the Monte Carlo code to take 70% of the systems assigned 
a period between 7 and 14 days and reassign them assuming the same power-law describing 
the rest of the distribution to periods between 3 and 5 days. The 70% reassignment fraction 
reflects a 90% completeness level and detection of one out of a probable 3 or 4 systems in this 
range. This approach approximates what would happen if binary systems with periods 7 and 
14 days migrated to shorter periods. The best fit using this approach is shown in Figure [13] 
and comes from reassigning systems with a period between 7 and 13.5 days to a period 
between 4 and 5 days (corresponding to the peak of the excess in the Cyg 0B2 distribution) 
with only a slight increase in /3 (0.5 and within the uncertainty of our results), a, 7, and 
the binary fraction are same as in Section 4. The normalized cumulative distribution for the 
Cyg 0B2 data is represented by the diamonds and the solid curve represents the normalized 
cumulative period distribution for the best fit using this version of the MC code. Though 
this fit requires a slightly more positive slope for the initial period distribution, the overall fit 
to the data is visually better than in Figure HJ and the multiple changes in the slope are also 
replicated well. A two-sided K-S test yields a compatibility of greater than 99% (compared 
to ~60% with the simple power- law approach). 

A plausible explanatio n for s uch period migration could be the mechanism proposed 



by iKrumholz fc Thompson! (120071 ) whereby the circularization of the orbit occurs primar- 
ily during the formation process instead. In Figure [H] we show the evolution of orbital 
period/semi-major axis and mass ratio for a binary system with a total mass of 30 M© as 
the primary loses mass to the secondary. The plotted curves show tracks of constant angular 
momentum as systems evolve toward shorter periods with mass ratios approaching unity, 
assuming that mass and angular momentum are conserved, 

a fqo\ f l + q 



«o V <? / V 1 + "Jo 



e.g., equation 8 of lKrumholz fc Thompson! 120071 ). Here Oq and go are the initial semi-major 



axis and mass ratio while a and q are the final semi- major axis and mass ratio after mass 
transfer. The two solid curves are tracks for systems with initial mass ratios g = 0.5 and 
initial periods of 7 and 14 days. These systems may end up having mass ratios near unity, 
but their final periods are longer than 5 days. The dotted curves are tracks for systems 
that begin with mass ratios q = 0.3 and periods of 7 and 14 days. These systems may 
end up with mass ratios near unity and periods in the range 3-5 days. This is precisely 
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where we observe an excess of systems with short periods in Cyg 0B2, and we also find 
a relative abundance of systems with mass ratios q =0.8-0.9 — not exactly "twins", but 
sys tems where significant rnass tra nsfer may have taken place. This figure indicates that, if 



the iKrumholz fc ThompsonI (120071 ) mechanism is the operative one producing an abundance 
of short-period, high-g systems, then some fraction of the population of 3-5 day systems 
should originate as systems with q ~ 0.3 with 7-14 day periods. While the evidence for 
"missing" 7-14 day systems in Cyg 0B2 is compelling, we cannot claim, on the basis of the 
present data, to see the corresponding lack of systems with q ~ 0.3. Our survey is sufficiently 
incomplete in this mass ratio range that observational biases preclude any firm conclusion. 



5.2. Role of Close Binaries in Energetic Phenomena 

The Cyg 0B2 radial velocity survey originated as an observational effort to inform 
theoretical models predicting the frequency of energetic events like supernova and 7-ray 
bursts that may have close binaries as progenitors. A key ingredient of these modeling 
efforts is the fraction of massive stars that might interact through mass transfer or undergo 
a common envelope phase, especially during the post-main-sequence evolution of either star 
in the system. R ed supergiants are the descendants of massive stars, so we adopt the RSG 



radii deduced by iLevesque et al.l ( 120051 ) which range between 200 R0 and 1500 R©. Their 
Figure 3 shows stellar evolutionary tracks indicating that stars in the mass range 15 - 30 
Mq may reach or exceed 1000 R© while stars with initial masses over 10 M0 can generally 
be expected to reach 500 R©. We adopt 500 R© as a maximum evolved radius for stars 
in our sample, meaning that systems with semi-major axes smaller than this value will 
interact strongly, probably through the formation of a com mon envelope. O ther systems 



may interact via mass transfer after Roche Lobe overflow. lEggletoru ( 1l983l ) provides an 
analytic description of the Roche Lobe radius as a function of mass ratio. His tabulated 
values show that the Roche Lobe radii for binary systems with mass ratios q =0.2 - 1.0 fall 
in the range 0.38a - 0.55a where a is the separation of the components. Consequently, we 
adopt 1000 Rq (4.7 A.U.) as a conservative threshold separation, below which systems are 
deemed to interact during the evolution of the system. In our Monte Carlo code, we flag 
all systems with periastron distances less than 1000 Rq as capable of interaction. Assuming 
a = 0.1, /3 = 0.2, 7 = —0.6, and Pmax = 10^ years, this predicts that 45% of the binary 
systems will have sufficiently close companions to interact during the binary's lifetime (37% 
if we adopt the slightly larger /3 = 0.5 from Section 5.1). 
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5.3. The Luminosity Contribution from Secondaries 

The luminous contribution from secondary stars, which are unresolved in photometric 
surveys, even with the HubbleSpaceTelescope or adaptive optics, is not negligible. Con- 
sequently, measuring accurate distances to star clusters using spectroscopic parallax, main 
sequence fitting, or similar standard-candle techniques is systematically in error if close com- 
panions are neglected. As part of our Monte Carlo code we tallied the mean V-band light con- 
tribution from primary and secondary stars by using a zero-age main sequence (ZAMS) rela- 



tion b etween stellar mass and absolute V magnitude (JMartins et al.ll2005 



Drilling fc Landolt 



2000). For the best fitting power law indices described above, allowing for periods as long 
as 10^ yr, we find that the secondary stars contribute 16% of the total stellar light at V- 
band. This is equivalent to a systematic error in the distance modulus of 0.16 mag and a 
systematic distance error of about 8%, in the sense that distances derived without correcting 
for luminous secondaries will be too small. We expect the secondary luminosity fractions 
in evolved populations to vary stochastically about this 16% ZAMS value as the primary 
and secondary components throughout the population become giants and supergiants. The 
intricacies of mass exchange between components will further limit any simple description 
of the time evolution in the nominal ZAMS value. 



6. Conclusions 

In this, our fifth work of this series, we present the results of a Monte Carlo analysis 
using a decade of radial velocity data on 114 massive stars in the Cyg 0B2 Association. 
This survey is unique in terms of number of massive stars and binary solutions ana l yzed. 



It irnproves on KF07 and other works (e.g., lEvans et al.l l2006l . ICarmany et al.lll980l . iGies 



19871 . iPourbaix et al.l |2004J ) in that we utilize the radial velocity information obtained over 
a larger time span and make use of the nearly complete orbital information — including 
eccentricities — on 22 of the 24 known massive binaries. The Cyg 0B2 binaries have been 
observed by spectroscopic and/or photometeric means and cover a period range of 1.47 days 
- 7.58 years (with all but a few having P < 26 days), an eccentricity range of 0-0.53, and a 
mass ratio range of 0.03-0.99. On the strength of these data, we compute the binary fraction 
for a large, young population (2-3 Myr) of massive stars from a common origin. In addition, 
we use this information to model the most probable intrinsic distributions for period, mass 
ratio, and eccentricity for short-period systems having P<26 days. 

• The data imply a hard lower limit on the binary fraction of 21% (i.e., 24/114 systems). 
The Monte Carlo analysis of the unbiased sample of 110 OB stars (i.e., not including 
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A36, A45, B17, and Schulte 5), indicates a massive star binary fraction of 44 ± 8%, 
assuming that the allowed period range is 1-1000 days (separations ~0.01 - few x 
10 AU). However, if we extrapolate the best fit power-law distributions over several or- 
ders of magnitude and allow for periods as long as 10^ years, the binary fraction may be 
as large as 90 ±10%. These results are consistent with the prob able number of binaries 



detect ed in the adaptive optics and HST survey of Cyg 0B2 by lCaballero-Nieves et al. 
( I2OIII ) and with the previous results of KF07. 



The best-fitting distributions of mass ratio, period, and eccentricity are described by 
power-laws with indices a = 0.1 ± 0.5, [3 = 0.2 ± 0.4, and 7 = —0.6 ± 0.3, respec- 
tively, assuming lower and upper bounds on the power-law distributions of q=0.005- 
1.0, P=l-1000 days, and e=0. 0001-0. 9. These values for a and (3 a re broadly consistent 



with t h e previous resul t s of K F07, as well as ISana et al.l (120081) . iKouwenhoven et al. 
( 20071 ). iGarmany et al.l ( 19801 ). and the analysis of ISana &: EvansI ( 201l[ ). The mass 



ratio distribution is no t cons istent with the "twin-heavy" population proposed by 
Pinsonneault fc Stanekl (|2006[ ). but suggests a weak preference for larger mass ratios 
(i.e., q > 0.8). 

The distribution of orbital periods is not consistent with any single power-law, owing 
to a significant excess of systems with periods around 3-5 days (peaking at 4-5 days) 
and a deficit of systems with periods around 7- 14 days. Such a 3-5 day excess has 
been observed among massive stars in NGC 6231 ( ISana et al.ll2008l ). but the latter phe- 
nomenon has not been previously reported. We propose that it may stem from binary 
systems with periods initially in the range of ~7-14 days migrating to shorter periods 
during a period of pre-main-se quence mass transfer, via the mechanism proposed by 
Krumholz fc Thompson! ( 120071 ). 



The fraction of massive-star binaries (or higher order systems) having companions 
close enough to interact during pre- or post-main-sequence evolution, defined here 
as requiring a semi-major axis <4.7 AU (~1000 R©), is ~45%. This fraction is a key 
ingredient in recipes that predict the cosmic rates of binary-induced phenomena arising 
from massive stars, such as some classes of supernova and 7-ray bursts. 



We present these findings in the anticipatio n that they, along with t he analyses of 
other clusters (s uch as the wo rk with Sco 0B2 bv IKouwenhoven et al.l 120071. NGC 6231 by 
Sana et al.ll20o3 . NGC 2244 bv lMahv et allbood . and NGC 6611 bv lSana et^l2009[ ). will aid 
our understanding of numerous issues involving massive stars, including their formation, the 
predicted rates of supernovae type Ib/c arising through binary channels, and whether there 
is a correlation between cluster density and the binary fraction therein. This ongoing survey 
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will continue to uncover massive spectroscopic binaries in Cyg 0B2. These new SBls and 
SB2s will help refine our characterization of the orbital parameter distributions and increase 
our understanding of the nature of the excess in orbital periods between ~3-5 days. While the 
present data has provided solid constraints on the massive binary characteristics in Cyg 0B2, 
it is possible that the binary fraction or distribution of periods/mass ratios/eccentricities 
varies temporally as a cluster ages or is a function of initial cluster density. Discerning 
whether these additional variables are important will require more extensive surveys over a 
range of age and formation environment. 
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from our referee in the process of making this a much better manuscript. We thank Mark 
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very appreciative for the continued advice and encouragement throughout the project from 
Ginny McSwain. Additionally, we are grateful for the continued support from the National 
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Table 1. Epheinerides for Cyg 0B2 Stars Used in Monte Carlo Analysis 



Date Vr (TV 

Star HJD-2,400,000 (km s^^) (km s^^) 



A36 


54,403.61 


-163.4 


5.6 


A36 


54,403.70 


-174.9 


5.0 


A36 


54,405.71 


103.4 


5.3 


A36 


54,406.63 


105.4 


3.7 


A36 


54,406.76 


108.6 


3.7 



N ote. — Column 1 gi ves the star designatio n in ei- 
therJSchuitc ( 1958) fSVComcr on et all J2OO2I ') (A,B) 
or iMasscv fc Thompson ( 1991) (MT) format. Col- 
umn 2 gives the heliocentric Julian date in the format 
HJD-2,400,000. Column 3 gives the cross-correlation 
velocity in km s~^. Column 4 gives the la cross- 
correlation error in km s~^. Table 1 is published in 
its entirety in the electronic edition of the Astrophys- 
ical Journal. A portion is shown here for guidance 
regarding its form and content. 
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Table 2. Massive Star Binary Statistics 



Cluster 


O Stars in 


Early-B Stars 


b.f. 


b.f. 


Reference 


ID 


Survey 


in Survey 


(O-Stars) 


(Massive Stars) 


ID 


NGC 2004 (LMC) 


4(1) 


83 (21) 


25% 


25% 


1 


NGC 330 (SMC) 


6(0) 


76(3) 


0% 


4% 


1 


NGC 346 (SMC) 


19(4) 


59 (19) 


21% 


29% 


1 


Trumpler 14 


7(1) 




14% 




2,3 


IC 1805 


10(2) 




20% 




4,5 


NGC 2244 


6 (1-2) 




17-33% 




6 


NGC 6231 


16 (10) 


33 (16) 


63% 


53% 


2,7 


NGC 6611 


9(4) 




44-67%'' 




2,3,8 


IC 2944 


16(7) 




44% 




2,3 


Trumpler 16 


20(7) 




35% 




2,3 


Cr 228 


21(5) 




24% 




2,3 


Cyg OB2 


44 (12) 


69 (8) 


27% 


18% 


9,10,11 



Note. — A selection of OB cluster/association studies examining the binary fraction of massive 
stars. The number of binaries are shown in parentheses. 



References. — flllEvans et al.l ||2006|): f2llGarciafc MermilliodI ||200J); (yiMermilliodI lll995ll: (4 ) 
iDe Becker et al.lll2006l '): f S'llRauw fc De Becked ||2004|'): (OllMahv et al.lll2009h i fT'tlSana et al.lll2008h : 
(8'l lSana et al.lll2009l) : ('9) lKiminki et al.l ||2008[) : (lOl lKiminki et al.lll2009h ; flD lKiminki et al.lll2010l l 

''Upper bound obtained from modeling complete fraction. 
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Fig. 1. — Histograms of the Cyg 0B2 orbital period distribution in log-Penod (days) (top 
panel), orbital eccentricity distribution (middle panel), and mas s ratio distributi o n (bo ttom 
panel). The orbital period and mass ratio distributions from iGarmany et al.l (jl980[ ) are 
illustrated by the dashed lines. 



-29- 



1.0 
0.8 

■f 0.6 

£ 0.4 

o 

o 

Ld 

0.2 

0.0 

1.0 

0.8 
_o 

2 0.6 

C/l 

in 

I 0.4 

0.2 

0.0 



o 






e 
h 



^ 



* 







o 



,o 



() 



® 



o. 



i \ 



O SB1 








• SB2 








)K9th 


SB 


Cata 


og 



10 15 20 

Period (days) 



25 



30 
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for 22 of the 24 known massive star binaries in Cyg 0B2. SB Is are represented by open 
circles and SB2s are represented by f illed circles. The dat a for massive star binaries of the 
Ninth Spectroscopic Binary Catalog (jPourbaix et al.ll2004J ) are over-plotted for comparison 
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Fig. 3. — Contour plots depicting the probabilities of K-S tests between the observed orbital 
parameter and radial velocity data of the "complete" sample and the MC code fits. Contours 
in the upper left panel are for a /3 = 0.2 and B.F. = 50%. Contours in the upper right panel 
are for a a = 0.1 and B.F. = 50%. Contours in the bottom left panel are for a 7 = —0.6 
and B.F. = 50%. Contours in the bottom right panel are for a a = 0.1 and 7 = —0.6. 
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Fig. 4. — Cumulative period distribution (diamonds) and conventional histogram for the 
observed orbital periods. The solid curve is the cumulative period distribution for the best 
MC code fit (/3 = 0.2). 
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Fig. 5. — Survey completeness as a function of orbital period for the best-fit power-law 
indices. 
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Fig. 6. — Normalized cumulative distribution (diamonds) and conventional histogram for 
the observed orbital eccentricities. The solid curve is the normalized cumulative distribution 
for the best MC code fit, 7 = -0.6. 
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Fig. 7. — Survey completeness as a function of eccentricity for the best-fit power-law indices. 
The solid curve represents the period range of 1-1000 days and the dot-dash curve represents 
the period range of 1-30 days. 
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Fig. 8. — Normalized cumulative distribution (diamonds) and conventional histogram for 
the observed orbital mass ratios. The solid curve is the normalized cumulative distribution 
for the best MC code fit, a = 0.1. 
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Fig. 9. — Survey completeness as a function of mass ratio for systems with periods P < 30 
days (dot-dash curve) and P < 1000 days (sohd curve) for the adopted best-fit power-law 
indices. . 
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Fig. 10. — x^ probability density function for radial velocity measurements from the Cyg 0B2 
survey. Systems with Pix^, v) <5% are probable velocity variables. 
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Fig. 11. — The cumulative distribution of x^ probabilities for the observed data (diamonds) 
and the best MC code fit (solid curve). The best fit represents a binary fraction of 44%. 
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Fig. 12. — Contour plots depicting the probabilities of K-S tests as in Figure El but for the 
observed orbital parameter and radial velocity data of the "unbiased" sample and the MC 
code fits and with B.F. = 44%. 



-40- 



c 
o 



J3 



Q 

(D 
_> 

-^^ 
O 

E 

D 

o 

CD 

"5 



1.0 



0.8 



0.6 



0.4 



0.2 



0.0 I 
0.0 













/o 








/ 


6 






^^ - 








^ — """"""^^ 






^ 


/ 


4 








/ 


f{log P) °< (Log P) 0-5 










i 


1 1 


2 




o 


/ 














r 


/ 


















^^ 


















o 


^ 


^ 


















E 

D 

cn 
o 

X 



0.2 0.4 0.6 0.8 1.0 

log Period (days) 



1.2 



1.4 



Fig. 13. — Normalized cumulative distribution (diamonds) and conventional histogram for 
the observed orbital periods. The solid curve is the cumulative distribution for the best MC 
code fit, assuming a 70% orbital period migration from P=7-13.5 days to P=4-5 days. 
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Fig. 14. — The evolution of orbital period/semi- major axis and mass ratio for a binary 
system with a total mass of 30 M© as the primary loses mass to the secondary. Arrows show 
the direction of evolution as systems shift toward shorter periods (smaller orbits) and larger 
mass ratio. The solid and dotted tracks show the evolution of systems with different initial 
masses ratios and orbital periods. Systems with initial periods P=7-14 days and q <0.3 end 
up as systems near 3-5 days and g ~1. 



